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Capillary-induced  transport  of  liquid  water  inside  the  porous  diffusion  media  (DM)  of  polymer  electrolyte 
fuel  cells  (PEFCs)  is  strongly  dependent  on  DM  pore  structure  and  material  properties.  As  such,  excessive 
liquid  in  the  DM  can  be  expelled  more  efficiently  into  flow  channels  by  proper  design  of  the  DM  structure. 
The  present  study  is  devoted  to  exploring  multiphase  transport  characteristics  by  considering  the  effects  of 
DM  pore  structure  and  material  properties.  Two  main  effects  on  overall  water  removability  are  examined, 
namely:  (i)  the  effect  of  immobile  liquid  saturation,  which  is  a  threshold  value  for  initiating  macroscopic 
capillary  transport  via  connected  small  liquid  droplets,  and  (ii)  the  effect  of  hydrophobic  spatial  variation, 
which  is  encountered  in  typical  DM  treated  with  PTFE.  Although  these  two  effects  are  expected  to  influence 
significantly  the  transport  characteristics  in  the  fuel  cell  DM,  they  have  been  ignored  in  most  two-phase 
fuel  cell  models  reported  in  the  literature.  In  the  present  work,  these  features  are  implemented  into  a 
one-dimensional,  multiphase  mixture  (M2 )  fuel  cell  model  along  the  through-plane  direction,  where  both 
anode  and  cathode  sides  and  the  membrane  are  fully  incorporated.  The  results  of  the  model  simulation 
clearly  demonstrate  the  dramatic  influence  of  the  amount  of  liquid  accumulation  and  capillary  transport 
characteristics  inside  the  DM.  The  findings  are  useful  for  designing  and  optimizing  DM  for  the  purpose  of 
effective  water  removal. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Water  management  in  polymer  electrolyte  fuel  cells  (PEFCs) 
affects  performance,  durability,  and  cold-start  characteristics.  Dif¬ 
fusion  media  (DM)  in  PEFCs  should  be  able  to  expel  efficiently 
excessive  water  condensed  in  the  catalyst  layer  (CL)  into  the  flow 
channels  in  order  to  minimize  flooding.  Typical  carbon-fibre-based 
DM  are  a  carbon-fiber-based  porous  materials  treated  with  PTFE, 
and  display  complex  pore  structures  and  non-uniform  wetting 
characteristics  [1-3].  Accordingly,  DM  pore  structure  and  material 
properties  are  key  factors  for  the  control  of  liquid  water  transport  in 
the  DM  and  thus  need  to  be  optimized  for  better  water  management 
of  PEFCs. 

Two-phase  transport  and  resulting  flooding  phenomena  in  DM 
have  been  investigated  using  both  modelling  and  experimental 
studies.  A  large  number  of  macroscopic  PEFC  models  reported  in 
the  literature  are  based  on  the  two-phase  Darcy’s  law  under  for  a 
homogeneous  pore  structure  and  uniform  wetting  characteristics 


*  Tel.:  +1  82  32  860  7312;  fax:  +1  82  32  868  1716. 
E-mail  address:  hcju@inha.ac.kr. 


[4-22].  Some  publications  further  analyze  the  effects  of  a  micro 
porous  layer  (MPL),  which  is  finer  and  more  hydrophobic  than  tra¬ 
ditional  DM  [18-21].  Although  macroscopic  fuel  cell  models  are 
based  on  descriptions  of  continuous  flow  and  transport,  the  effect 
of  immobile  liquid  saturation,  which  constitutes  a  threshold  for 
a  continuous  liquid  phase,  has  been  largely  ignored.  Only  a  few 
researchers  [20-22]  have  examined  its  influence  on  the  amount  of 
liquid  accumulation  in  DM. 

Simultaneously,  efforts  have  also  been  made  to  investigate 
experimentally  the  effects  of  DM  pore  structure  and  hydrophobic- 
ity  on  water  transport  through  DM  as  well  as  on  cell  performance. 
Kong  et  al.  [23]  showed  the  influence  of  DM  pore-size  distribution 
on  cell  I-V  performance  curves.  Park  et  al.  [24]  and  Qi  and  Kauf¬ 
man  [25]  compared  and  analyzed  cell  performance  using  DM  with 
different  pore  structures,  PTFE  contents,  and  MPLs.  In  addition,  a 
neutron  radiography  technique  was  employed  to  examine  the  effect 
of  DM  properties  on  liquid  accumulation  in  PEFCs  [26-28].  These 
experimental  studies  provide  strong  evidence  that  DM  structure 
and  wettability  affect  water  transport  characteristics  in  fuel  cell 
DM. 

Macroscopic  models  have  inherent  limitations  in  accurately 
describing  multi-phase  flows  through  heterogeneous  porous 
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Nomenclature 

a  water  activity  or  effective  catalyst  area  per  unit  of 

total  volume  (m2  m-3) 

A  area(m2) 

Dk  mass  diffusivity  of  species  k  (m2  s-1 ) 

EW  equivalent  weight  of  dry  membrane  (kg/mol) 

F  Faraday  constant  (96,487  C  mol-1 ) 

I  current  density  (A  m-2 ) 

jl  diffusive  mass  flux  of  i  phase  (kg  m-2  s-1 ) 

J  Leverett  function 

kr  relative  permeability 

I<  hydraulic  permeability  ( m2 ) 

m  mass  fraction 

M  molecular  weight  (kg  mol-1 ) 

n  number  of  electrons  in  electrochemical  reaction  or 

diffusivity  correction  factor 
nd  electroosmotic  drag  coefficient 

P  pressure  (Pa) 

Pc  capillary  pressure  (Pa) 

RH  inlet  relative  humidification 

Ru  universal  gas  constant  (8.314J  mol-1  K-1 ) 

s  liquid  saturation 

sim  immobile  liquid  saturation 

sint  interfacial  liquid  coverage 

sr  reduced  liquid  saturation 

T  temperature  (K) 

u  fluid  velocity  and  superficial  velocity  in  porous 

medium  (ms-1) 

V  volume  (m3) 

Greek  letters 

a  transfer  coefficient 

y  advection  correction  factor 

8i  thickness  of  component  i 

s  volume  fraction  of  gaseous  phase  in  porous  region 

0  contact  angle  (°) 

A  membrane  water  content  ( mol  H2  O/mol  S03  - ) 

Aa  relative  mobility  of  phase  a 

fi  viscosity  (kg  m-1  s-1 ) 

p  density  (kg  m-3) 

pmem  dry  membrane  density  ( kg  m-3 ) 
a  surface  tension  (Nm-1)  or  electronic  conductivity 

(Sm-1) 

r  viscous  shear  stress  (N  m-2) 

v  kinematic  viscosity  (m2  s-1 ) 

Subscripts 
a  anode 

avg  average  value 

c  cathode 

CL  catalyst  layer 

DM  diffusion  medium 

g  gas  phase 

GC  gas  channel 

H2  hydrogen 

i  species  index 

mem  membrane 

N2  nitrogen 

02  oxygen 

sat  saturation  value 

w  water 

0  standard  condition,  298.15  K  and  101.3  kPa  (1  atm.) 


Superscripts 

eff  effective  value  in  porous  region 

mem  membrane 

g  gas 

1  liquid 

sat  saturation  value 


media.  In  response  to  this,  pore-scale  modelling  and  simulations 
have  recently  been  conducted  by  several  fuel  cell  research  groups, 
with  the  aim  of  more  precisely  investigating  the  liquid  water 
transport  mechanism  through  a  fibrous  DM  structure  [29-34].  In 
addition  to  providing  constitutive  relations  under  a  carbon-fibre- 
based  DM  microstructure  (e.g.,  capillary  pressure  and  relative 
permeability  as  a  function  of  liquid  saturation),  which  are  employed 
as  inputs  for  macroscopic  two-phase  fuel  cell  models,  the  pore- 
level  simulations  revealed  that  pore  morphology  and  non-uniform 
spatial  distribution  of  wettability  strongly  influence  the  capillary 
transport  characteristics  and  liquid  water  accumulation  inside  DM. 
While  pore-level  models  are  capable  of  predicting  more  realis¬ 
tic  liquid  water  profiles,  these  models  are  computationally  more 
demanding  compared  with  macroscopic  models  based  on  the  vol¬ 
ume  averaging  theory.  Furthermore,  they  are  more  difficult  to 
incorporate  into  a  comprehensive  fuel  cell  model. 

In  the  present  work,  we  describe  the  effects  of  immobile  liq¬ 
uid  saturation  and  non-uniform  DM  wettability  using  macroscopic 
two-phase  fuel  cell  simulations.  These  features  are  implemented 
into  a  two-phase  PEFC  model.  The  model  is  a  one-dimensional, 
multiphase  mixture  (M2)  fuel  cell  model  along  the  through-plane 
direction  that  fully  incorporates  both  the  anode  and  cathode  sides 
and  the  membrane.  Model  simulations  and  analyses  yield  a  better 
macroscopic  understanding  of  two-phase  transport  processes  and 
flooding  physics.  This  provides  useful  information  in  designing  and 
optimizing  DM  for  the  purpose  of  effective  water  removal. 


2.  Two-phase  PEFC  model 

2.1.  Multi-phase  mixture  (M2 )  model 

The  one-dimensional,  two-phase  flow  model  presented  in  the 
current  study  is  based  on  a  percolation  theory  proposed  by 
Pasaogullari  and  Wang  [4]  and  Nam  and  Kaviany  [20].  In  the  the¬ 
ory,  liquid  water  droplets  that  preferentially  occupy  larger  pores  in 
hydrophobic  porous  DM  can  agglomerate  and  become  a  continu¬ 
ous  liquid  phase.  The  continuous  phase  enables  the  liquid  water 
to  move  through  the  DM  under  a  liquid  pressure  gradient.  The 
macroscopic  multi-phase  flow  phenomena  through  porous  media 
are  traditionally  modelled  by  a  multiphase  mixture  (M2)  approach 
where  consideration  of  immobile  liquid  saturation  is  indispensable 
for  the  description  of  continuous  flow  and  transport.  In  the  present 
study,  a  previously  described  M2- based  fuel  cell  model  [4,5,12,18] 
is  reformulated  by  adding  a  new  feature,  namely,  immobile  liquid 
saturation. 

In  order  to  quantify  individual  phase  distributions  in  porous 
media,  liquid  saturation,  s,  is  defined  as  the  volume  fraction  of  pores 
occupied  by  the  liquid  phase: 


s  = 


V1 

V 


(1) 


When  the  immobile  liquid  saturation,  sim,  is  taken  into  account, 
a  reduced  liquid  saturation,  sr,  which  can  be  interpreted  as  the 
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fraction  of  transportable  liquid,  has  to  be  introduced  as  follows: 


sr  = 


s-% 

1  Sir 


(2) 


Since  two-phase  transport  inside  porous  DM  is  mainly  governed  by 
capillary  motion  between  liquid  and  gas  phases,  capillary  pressure, 
defined  as  the  difference  between  gas  and  liquid  pressures,  is  a  basic 
parameter  in  the  M2  model.  The  capillary  pressure,  Pc,  is  expressed 
as  [35]: 


PC  =  PS  -Pl  =GCOS0{  - 


(f) 


1/2 


J(ST 


(3) 


where  e,  K ,  and  9  are  the  porosity,  the  permeability,  and  the  contact 
angle  of  the  given  porous  media,  respectively.  In  Eq.  (3),  the  Leverett 
function,  J(sr),  meaning  a  dimensionless  capillary  pressure,  can  be 
expressed  as  a  function  of  the  reduced  liquid  saturation,  sr,  for  the 
given  hydrophobic  DM: 


J{st 


if  s  <  sim 


417sr-  2.120s?  +  1.263s?  if  S  >  S 


(4) 


It  should  be  noted  here  that  the  standard  Leverett  J-function  form 
in  Eq.  (4)  was  formulated  to  characterize  the  liquid  transport 
in  more  uniform  porous  media  such  as  soils  and  rocks  [36,37]. 
A  more  realistic  capillary  pressure-saturation  study  for  real  fuel 
cell  DM  was  recently  conducted  by  Gostick  et  al.  [3]  and  Kum- 
bur  et  al.  [2,39,38].  In  particular,  Kumbur  et  al.  [39]  conducted 
capillary  pressure-saturation  measurements  under  typical  fuel¬ 
cell  compressed  DM  and  operating  temperatures  (20,  50,  and 
80  °C).  Experimental  measurements  were  compared  with  capil¬ 
lary  pressure-saturation  correlation  test  data  obtained  with  the 
traditional  Leverett  function.  According  to  the  comparison,  the 
deviation  between  the  Leverett  function  and  new  correlations  was 
not  significant,  especially  for  the  practical  range  of  liquid  satura¬ 
tion  (0  <  s  <  0.5)  in  DM  during  fuel  cell  operation.  Therefore,  we  have 
employed  the  standard  Leverett J-function,  i.e.,  Eq.  (4),  in  this  study. 

On  the  other  hand,  a  relative  permeability  term  should  be  intro¬ 
duced  for  the  two-phase  porous  region  because  the  available  pore 
space  in  porous  DM  is  shared  by  both  gas  and  liquid  phases.  The 
relative  permeability  for  an  individual  phase  can  also  be  expressed 
in  terms  of  reduced  liquid  saturation,  sr,  as  follows: 


kl  = 


0 


if  0  <  S  <  Sir 


=  (S-Sim  y 

U -W 


if  S  >  Sir 


k?  =  (l-s)3 


(5a) 


(5b) 


where  k\  and  kf  denote  the  relative  permeabilities  of  non-wetting 
liquid  and  wetting  gas  phases,  respectively,  for  hydrophobic  DM. 
Fig.  1  schematically  shows  the  relative  permeability  curves  as  a 
function  of  liquid  saturation,  s,  as  defined  in  Eq.  (1).  It  is  seen  that 
the  relative  permeability  of  the  non-wetting  liquid  phase,  k\  is  zero 
until  liquid  saturation,  s  exceeds  the  threshold  value  at  immobile 
liquid  saturation. 

Finally,  all  mixture  properties  in  the  M2  formulation  are  defined 
as  follows  [35]. 

Mixture  density: 


p  =  pxs  +  pg{  1  -s) 


(6) 


where  the  gas  mixture  density,  pg,  described  by  the  ideal  gas  law 
varies  with  its  composition  (mass  fractions,  mf).  That  is, 


"-(sr) 


1 


E  i™VMi 


(7) 


0.0 


1.0 


Liquid  saturaiton,  s 


Fig.  1.  Relative  permeabilities  for  non-wetting  liquid  and  wetting  gas  phases  as 
function  of  liquid  saturation  for  hydrophobic  DM.  Dashed  line  denotes  immobile 
liquid  saturation. 


Mixture  velocity: 
pxux  +  pgug 


u  = 

P 

Species  mass  fraction: 

pxsmx  +  pg(  1  -  s)mf 

m,  = - - - — 

P 

Kinematic  viscosity: 

(k\  k?N_1 


(8) 


(8) 


(9) 


where  vg  is  the  kinematic  viscosity  of  a  gas  mixture  varying  with 
gas  composition  [40]: 


.  _  1  y-  xum 

pS  pS  Z-/  ’ 


where  0,,  =  — — 

J  V8 


1  + 


-1/2 


1  + 


\  1/2 

P±  \ 

HJ 


[Mi 


1/4' 


(10) 


and 


fii  (Ns/m2)  =  < 


/zh2  =  0.21  x  io-6T0-66 
/xw  —  0.00584  x  io_6T1,29 
/xN2  =  0.237  x  l0-6r°-76 
/Xq2  =  0.246  x  IQ-6!0  78 


T  in  K 


Relative  mobility: 

I'.Sw 

V1 

Xs  =  1  -  A1 

2.2.  Fuel  cell  model  description  and  assumptions 


(11) 

(12) 


In  this  study,  for  complicity,  the  computational  domain  is  con¬ 
fined  to  the  anode  DM,  the  catalyst  coated  membrane  (CCM)  and  the 
cathode  DM.  Therefore,  both  the  anode  and  cathode  gas  channels 
(GCs)  are  excluded  from  the  computational  domain  by  applying 
proper  boundary  conditions  to  the  DM|GC  interfaces.  In  addition, 
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Anode 
GC 
Electro- Os 

Membrane 
or  water  coi 


H2  consumptio 


Cathode 

GC 


Water  production  by  ORR  : 


2  F 


02  consumptio  i  by  ORR  : 


x  (Thru-plane) 


Fig.  2.  Schematic  of  modelling  domain  and  transport  phenomena,  where  HOR  and  ORR  denote  hydrogen  oxidation  reaction  in  anode  side  and  oxygen  reduction  reaction  in 
the  cathode  side,  respectively. 


the  cell  is  assumed  to  be  under  steady-state  and  isothermal  condi¬ 
tions  in  order  to  simplify  numerical  analysis  and  emphasize  liquid 
water  formation.  The  schematic  diagram  presented  in  Fig.  2  depicts 
the  associated  transport  processes  in  this  model. 


2.3.  Conservation  equations 


Mass  conservation  in  steady-state  along  the  through-plane 
direction  (x)  is 

^(pu)  =  0  (13) 

If  Eq.  (13)  is  integrated  along  the  x  direction,  the  total  mass  flux 
through  the  porous  DM  is  given  by  the  following  relationship. 
Mass  transport  in  the  anode  DM: 


Pu  =  ^mh2  + 


j  nmem  a  \  > 

n  _ _ nmem  P _ 

dF  w  EW  dx. 


Mw  - 


Kmem  dpl 

v1  dx 


mem 

(14) 


Mass  transport  in  the  cathode  DM: 


The  equation  can  be  written  in  one-dimensional  form  as 


JCdP 

pu=  -  — 

v  dx 


(16) 


The  steady-state,  one-dimensional  species  conservation  equation 
can  be  written  as  [35]: 


_d 

dx 


( K;pm,u ) 


d_ 

dx 


m\)f] 


(17) 


where  jl  in  the  second  term  on  the  right-hand  side  of  Eq.  (17) 
denotes  the  diffusive  mass  flux  of  liquid  phase  relative  to  the  whole 
multiphase  mixture  [35].  It  can  be  expressed  in  a  one-dimensional 
form  as  follows: 


j1  =  plu[  -  Xxpu 


(18) 


On  the  other  hand,  the  term  on  the  left-hand  side  of  Eq.  (17)  repre¬ 
sents  the  advective  term,  where  the  advection  correction  factor,  y* 
is  given  by  [35]: 

p{Xxmx  +  X%mg) 

Yi  =  7  x  \  ' - (19) 

(s/rmj  +  (1  -  s)pSmf ) 


pti  —  ——Mq2  +  —  Mw 


mem 

(15) 

where  the  terms  on  the  right-hand  side  of  Eqs.  (14)  and  (15)  rep¬ 
resent  the  reactants  and  product  fluxes  due  to  the  electrochemical 
reactions  and  water  flux  across  the  membrane  that  is  driven  by 
the  hydraulic  pressure  gradient,  diffusion  due  to  the  water  content 
gradient,  and  electro-osmotic  drag  due  to  proton  flux  (see  Fig.  2). 

For  the  momentum  conservation  equations  for  flow  through  the 
porous  DM,  Darcy’s  law  is  used  instead  Navier-Stokes  equations. 
For  the  two-phase  region,  the  Darcy’s  equations  are  derived  accord¬ 
ing  to  the  M2  model  for  both  single-  and  two-phase  porous  regimes. 


H)-l 


_  \  J _ nmemr _ 

F  w  EW 


dA2 

dx. 


Mw- 


7Cmem  dP1 
v1  dx 


Species  diffusivity  in  the  gas  mixture,  Df ,  in  the  first  term  on  the 
right-hand  side  of  Eq.  (17),  is  defined  such  that  summation  of  inter¬ 
species  diffusion  within  the  gas  phase  will  be  equal  to  zero  [40]: 


Dg  = 


1  -  Xj 


Zjxj/Di,. 


1/2 


u  „  1.013  x  10_7T175  /  1  1 

where  D{  /  = - =-  ttt  +  tt 

p(x,,/3+x;/3)  vMi  Mj 

XHW  =  7.07,  Xw  =  12.7,  Xnw  =  17-9,  Xow  =  l6-6 


(20) 


Note  that  gaseous  diffusive  transport  for  a  porous  medium  can  be 
controlled  by  the  Knudsen  diffusion  effect  due  to  molecule-to-wall 
collisions,  as  well  as  molecular  diffusion  caused  by  molecular  col- 
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lisions,  as  described  in  Eq.  (20).  The  Knudsen  diffusion  coefficient 
can  be  computed  according  to  the  kinetic  theory  of  gases  as  follows: 


_2/8 RuT\ 

3  V  nMi  )  p 


(21) 


The  effective  diffusivity  of  species,  i  in  the  gas  mixture,  is  then 
obtained  by  combining  both  molecular  and  Knudsen  diffusion 
effects  with  the  effects  of  porosity  and  tortuosity  of  the  porous 
medium  by  means  of  the  Bruggeman  correlation  [41  ]: 


D?’eff  =  [e(l  -  s)]n 


1  1 

D?  + 


-l 


(22) 


where  it  is  seen  that  the  effective  gas  diffusivity  Df,eff  in  the  two- 
phase  region  is  a  function  of  both  the  porosity,  £,  and  the  liquid 
saturation,  s. 

Eq.  (17)  can  be  written  for  water  species  as 


_d 

dx 


(ywpmwu)  — 


d_ 

dx 


ml!)'1] 


(23) 


Again,  this  equation  is  integrated  over  the  porous  DM  to  give: 


water  transport  in  the  anode  DM: 

ywpmwu  =  pgD^eff^(m®,)  +  (mg  -  m^)/1 


+ 


'  1  nmem  Pmem  dA  > 

d F  w  EW  dx; 


Kmem  dpl 

Mw - ; - -j— 

v1  dx 


water  transport  in  the  cathode  DM: 


mem 

(24) 


ywpmwu  =  /OgDg;eff^(mg,)  +  (mg  -  m'w)j 1 


+ 


/(mem  dP1 


n  j  nmemP^  dA.  > 
F  w  EW  dx  j 


Myv 


dx 


(25) 


As  noted  earlier,  this  study  takes  into  account  the  water  transport 
across  the  electrolyte  phase  driven  by  three  modes  [the  hydraulic 
pressure  gradient,  the  gradient  of  membrane  water  content,  and 
the  electro-osmotic  drag  due  to  proton  flux]  can  be  expressed  by 
the  following  equation: 
water  transport  across  membrane: 


d_ 

dx 


pmem  gm  dA  \  d 

EW  Dw  dx  J  Mw  dx 


Kmem  dpl  \ 

~v^~~dx  J 


=  0 
(26) 


The  transport  properties  of  electrolytes  are  correlated  with  the 
water  content  of  the  membrane,  A, which  is  in  turn  a  function  of 
the  water  activity,  a,  as  follows  [42]: 


a  = 


C*RUT 

Psat 


(27) 


f  Ag=0.043  +  17.81a  -  39.85a2  +  36.0a3  for  0  <  a  <  1 
\A'  =  22  (28) 


The  electro-osmotic  drag  coefficient,  nd,  and  the  water  diffu¬ 
sion  coefficient  in  the  membrane,  D™em  have  been  described  by 
Zawodzinski  et  al.  [43]  and  Motupally  et  al.  [44]: 


fl  A<A§  (a  =  l) 
d  |  2.5  A  =  A1 


(29) 


f  3.1  X  10“7A  (e°-28i  -  i)e(-2346/T)  for  0  <  A  <  3 
4.17  x  10_8A  (1  +  161  e_7-)e<_2346/r)  otherwise 

(30) 


Meanwhile,  a  one-dimensional  oxygen  species  conservation  equa¬ 
tion  provided  by  the  M2  formulation  is  obtained  from  Eq.  (17): 


d  ,  ,  d 

^yo2pm02u)=  ^ 


dx 


dx 


^S<f>So2) 


+  ^K2j‘] 


(31) 


In  this  study,  no  oxygen  is  assumed  in  the  anode  side,  and  hence 
Eq.  (31)  can  be  integrated  over  the  cathode  DM  to  give: 
oxygen  transport  in  cathode  DM: 


yo2pm02u  =  )  +  ms0J  ~  j^Mo. 


4  F 


(32) 


2.4.  Boundary  conditions 


The  aforementioned  one-dimensional  model  requires  boundary 
conditions  for  water  and  oxygen  species  at  the  GC|DM  inter¬ 
faces.  The  interfacial  mass  fractions  for  oxygen  and  water  species, 
mi ,dm/go  are  determined  by  the  cell  operating  pressure  (P),  temper¬ 
ature  (T),  GC  relative  humidity  (RHGc),  and  interfacial  liquid  droplet 
coverage  (sint)  due  to  channel  flooding  conditions: 

mw,DM/Gc(Sint>  RHgC?  T,  P) 

_  P^int  +  (1  —  5[nt)(PwMw/RuT) 

P^int  +  PS(^  —  ^int) 

_  P^int  +  (1  —  5int)((RHGcPsatMw)/PuT)  (33) 

pl  Sint  +  P^O  —  Sjnt) 


rn02,DM/Gc(sint^  RH,  T,  P) 

_  (1  “  sint)^02(^02/Ru T) 
plaint  +  Ps(l  —  Sint) 

_  (1  —  Sint)Mo2(P  —  Pw  —  Pn2/RuT)  ^4 j 

P^Sint  +  PS(^  —  Sjnt) 

2.5.  Numerical  procedures 

The  first-order  ordinary  differential  equations  derived  in  the 
foregoing  section  are  solved  separately  in  three  different  regions, 
i.e.,  anode  DM,  CCM,  and  cathode  DM.  In  order  to  connect  these 
differential  equations  for  the  three  domains,  an  initial  estimated 
value  is  provided  for  the  water  flux  across  the  membrane  to  set  up 
the  interfacial  boundary  conditions.  An  iterative  procedure  is  used 
to  improve  the  initial  estimate  where  the  individually  calculated 
interfacial  fluxes  in  adjoining  domains  should  be  matched  at  the 
interface.  The  iterations  proceed  until  the  relative  error  falls  below 
the  convergence  criterion  (10-7). 

3.  Results  and  discussion 

In  this  study,  the  aforementioned  one-dimensional  PEFC  model 
is  applied  to  typical  fuel  cell  geometry.  The  cell  component  proper¬ 
ties  and  physical  parameters  for  the  one-dimensional  simulations 
are  shown  in  Table  1.  It  is  assumed  for  all  simulation  cases  that  the 
cell  is  operated  at  a  constant  temperature  of  80  °C  and  a  pressure 
of  1.5  atm.,  which  are  typical  of  practical  PEFC  stack  operations.  In 
addition,  in  order  to  focus  on  the  cell  flooding  phenomena,  GC  rela¬ 
tive  humidity  (RHGG)  is  assumed  to  be  unity,  which  means  that  the 
gas  stream  in  the  flow  channel  is  fully  saturated. 

The  effect  of  immobile  liquid  saturation  on  the  DM  liquid  water 
distribution  is  illustrated  in  Fig.  3.  The  five  cases  shown  correspond 
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Table  1 

Material  properties  and  physical  parameters 


Description 

Value 

Porosity  of  DM,  £Dm 

0.6 

Permeability  of  DM,  7CDm 

1.0  x  10-12  m2 

Thickness  of  DM,  5dm 

0.3  mm 

Hydraulic  permeability  of  membrane,  I<mem 

5.0  x  10-20  m2 

Thickness  of  membrane,  5mem 

0.018  mm 

Dry  membrane  density,  pmem 

2000  kg  nrr3 

Equivalent  weight  of  electrolyte  in  membrane,  EW 

1.1  kg  mol-1 

Faraday  constant,  F 

96,487  C  mol-1 

Universal  gas  constant,  Ru 

8.314J  molK-1 

Surface  tension,  o 

0.0625  Nm-1 

Liquid  water  density,  p1  (80  °C) 

972  kg  irr3 

Liquid  water  viscosity,  p,1 

3.5  x  10-4  Ns m-2 

to  different  immobile  liquid  saturation  values  that  range  from  0  to 
0.2,  where  the  immobile  liquid  saturation  for  each  case  is  assumed 
to  be  constant  throughout  the  entire  DM,  although  it  is  likely  to 
vary  in  the  random  fibrous  DM.  These  cases  assume  no  interfacial 
liquid  coverage  (sic  =  0)  at  the  DM  surface.  This  means  that  most  of 
the  liquid  water  transported  through  the  DM  is  evaporated  at  the 
interface  of  the  DM  and  GC.  Thus,  channel  flooding  is  minimal  [12]. 
It  is  clearly  seen  in  Fig.  3  that  the  amount  of  liquid  water  accumula¬ 
tion  in  the  cathode  DM  is  strongly  influenced  by  the  magnitude  of 
immobile  liquid  saturation.  Without  the  effect  of  immobile  liquid 
saturation,  i.e.,  sim  =  0,  the  overall  liquid  saturation,  savg,  averaged 
along  the  DM  thickness  is  0.138  (13.8%)  and  the  maximum  liquid 
saturation  of  0.1 58  ( 1 5.8%)  is  predicted  to  be  near  the  DM  |  CCM  inter¬ 
face.  As  the  immobile  liquid  saturation  increases,  both  the  overall 
and  maximum  liquid  saturations  of  the  DM  increase.  When  the 
immobile  liquid  saturation  is  assumed  to  be  0.2,  the  overall  and 
maximum  liquid  saturations  are  predicted  to  be  0.295  and  0.320, 
respectively.  Generally,  the  immobile  liquid  saturation  of  porous 
media  is  increased  when  the  spatial  heterogeneity  in  porous  media 
also  increases.  Therefore,  it  can  be  concluded  from  Fig.  3  that  the 
level  of  DM  flooding  is  strongly  influenced  by  the  spatial  hetero¬ 
geneity  of  the  DM. 

Fig.  4  displays  the  effect  of  interfacial  liquid  coverage  on  cath¬ 
ode  DM  liquid  saturation  profiles,  where  sic  =  20%  in  Cases  2  and 


interface  interface 

Fig.  3.  Liquid  saturation  profiles  in  cathode  DM  as  function  of  immobile  liquid  sat¬ 
uration,  sim  (7=1.5  A  cm-2,  DM  contact  angle,  0  =  110°,  sic  =  0  for  all  cases).  Overall 
liquid  saturations,  savg,  i.e.,  averaged  along  DM  thickness  for  cases  of  Sim  =  0.0,  0.05, 
0.10,  0.15,  and  0.20  are  0.138,  0.166,  0.209,  0.252,  and  0.295,  respectively. 


interface  interface 

Fig.  4.  Liquid  saturation  profiles  in  cathode  DM  with  different  immobile  liquid  sat¬ 
uration,  sim,  interfacial  liquid  coverage,  sic  (7=  1.5  A  cm-2,  DM  contact  angle,  0  =  110° ). 
Case  1 :  sim  =  0%  and  sic  =  0% ;  case  2 :  sim  =  0%  and  sic  =  20% ;  case  3 :  sim  =  10%  and  sic  =  0%; 
case  4:  sim  =  10%  and  sic  =  20%. 

4  means  that  the  20%  of  the  cathode  DM  interfacial  surface  is 
covered  by  liquid  water  droplets.  This  implies  that  appreciable 
channel  flooding  exists  for  these  cases.  Comparing  cases  1  and  2 
which  have  no  immobile  liquid  saturation  effect,  Sjm  =  0%,  much 
higher  liquid  saturation  is  predicted  in  case  2  which  causes  a 
large  difference  in  average  liquid  saturation,  savg  (Asavg  =  6.9%). 
This  indicates  that  interfacial  coverage  has  a  significant  effect  if 
immobile  liquid  saturation  is  assumed  to  be  small.  On  the  other 
hand,  when  immobile  liquid  saturation  of  10%  is  considered,  the 
difference  in  the  average  liquid  saturation  between  cases  3  and  4 
is  smaller  (Asavg  =  1.6%).  The  comparison  of  the  four  curves  pre- 


DM/GC  Fractional  distance  DM/CCM 

interface  interface 

Fig.  5.  Liquid  saturation  profiles  in  cathode  DM  with  different  current  density, 
7,  immobile  liquid  saturation,  Sim,  and  interfacial  liquid  coverage,  SiC.  Case  1: 
7  =  1.5  A  cm-2 ,  sim  =  0%,  and  sic  =  0%;  case  2 :  7  =  0.8  A  cm-2 ,  sim  =  0%,  and  sic  =  10%;  case 
3:  7= 0.4  A  cm-2,  sim  =  0%,  and  sic  =  20%;  case  4:  7=  1.5  A  cm-2,  sim  =  10%,  and  sic  =  0%; 
case  5:  7  =  0.8Acnrr2,  sim  =  10%,  and  sic  =  10%;  case  6:  7  =  0.4Acnrr2,  sim  =  10%,  and 
sic  =  20%.  DM  contact  angle,  0  is  assumed  to  be  110°  for  all  cases. 
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CCM/DM  Distance  from  MEM,  [pm]  DM/GC 

interface  interface 


Fig.  6.  Liquid  saturation  profiles  in  cathode  DM  with  variation  of  DM  contact  angle. 

sented  in  Fig.  4  clearly  demonstrates  that  the  effect  of  immobile 
liquid  saturation  is  less  significant  with  a  high  level  of  DM  flood¬ 
ing,  e.g.,  the  presence  of  liquid  coverage  on  the  DM  interfacial 
surface. 

Fig.  5  gives  a  plot  of  liquid  saturation  profiles  in  the  cathode  DM 
at  three  different  current  densities,  0.4, 0.8,  and  1.5  A  cm-2,  where  a 
higher  value  of  interfacial  coverage  is  applied  to  the  lower  current 
density  cases.  Therefore,  the  assumed  values  of  interfacial  cover¬ 
age  at  0.4,  0.8,  and  1.5  A  cm-2  are,  respectively,  20, 10,  and  0%.  The 
assumption  is  based  on  experimental  observations  that  the  droplet 
detachment  diameter  from  the  DM  interfacial  surface  is  inversely 
proportional  to  the  air  flow  velocity  in  the  channels,  and  thus  to  the 
current  density  under  the  given  air  stoichiometry  [45].  Two  cases 
are  considered  for  each  current  density.  One  is  without  the  effect  of 
immobile  liquid  saturation  (cases  1  -3  for  sim  =  0.0)  and  the  other  is 
with  consideration  of  10%  immobile  liquid  saturation  (cases  4-6  for 
Sim  =  0.1 ).  From  a  comparison  of  these  cases,  it  is  clearly  seen  that 
the  effect  of  immobile  liquid  saturation  is  smaller  at  lower  current 
density.  The  smallest  effect  is  predicted  at  0.4  A  cm-2,  due  to  the 
larger  interfacial  coverage,  i.e.,  20%. 

It  is  expected  that  the  wetting  characteristics  of  DM  vary  spa¬ 
tially  due  to  heterogeneous  DM  pore  structures  as  well  as  anomalies 
in  the  PTFE  treatment  of  DM.  The  effects  of  spatial  variation  of  DM 
wettability  (contact  angle,  #)  on  the  liquid  water  distribution  in 
the  DM  are  investigated  in  Fig.  6.  According  to  Eq.  (3),  the  capil¬ 


lary  pressure  Pc  is  directly  proportional  to  the  cosine  of  the  DM 
contact  angle  (cos#).  Therefore,  the  linear  variation  of  cos#  along 
the  DM  thickness  is  considered,  as  shown  schematically  in  Fig.  6a. 
Case  1  is  a  baseline  case  without  variation  of  cos#  (uniform  con¬ 
tact  angle).  The  absolute  value  of  cos#  linearly  increases  (case 
2)  or  decreases  (case  3)  from  the  CCM|DM  interface  towards  GC, 
based  on  the  given  ratio  of  cos  #  in  Fig.  6a.  Liquid  saturation  dis¬ 
tributions  for  the  three  cases  are  shown  in  Fig.  6b.  It  is  clearly 
seen  that  the  direction  of  contact  angle  variation  along  the  DM 
thickness  strongly  affects  both  the  amount  of  liquid  water  accu¬ 
mulation  and  the  shape  of  the  liquid  saturation  profiles  inside  the 
DM.  Higher  liquid  saturation  is  predicted  in  case  2  where  the  DM 
hydrophobicity  increases  towards  GC.  This  indicates  that  liquid 
water  tends  to  be  more  easily  accumulated  in  relatively  hydrophilic 
pores  rather  than  being  transported  towards  relatively  hydrophobic 
pores.  On  the  other  hand,  case  3  has  lower  liquid  water  accumu¬ 
lation,  because  liquid  water  removal  inside  the  DM  is  facilitated 
when  the  DM  hydrophobicity  decreases  towards  GC.  The  liquid 
water  transport  characteristics  also  influence  the  shape  of  the  liq¬ 
uid  saturation  profile  in  the  DM,  where  the  concave-shaped  liquid 
saturation  profile  for  case  2  differs  from  the  convex  shape  for  cases 
1  and  3. 

4.  Conclusions 

A  one-dimensional,  multiphase  mixture  (M2)-based  fuel  cell 
model  has  been  developed.  The  model  takes  into  account  two 
key  features  in  order  to  capture  more  accurately  the  accumula¬ 
tion  and  distribution  of  liquid  water  inside  DM.  The  first  feature  is 
the  effect  of  immobile  liquid  saturation,  which  is  a  threshold  value 
for  continuous  liquid  phase,  and  the  second  feature  is  the  effect 
of  spatial  variation  of  DM  wettability,  which  is  mainly  caused  by 
the  heterogeneity  of  DM  pore  structures  as  well  as  non-uniform 
PTFE  treatment  of  DM.  One-dimensional  PEFC  simulations  have 
been  performed  to  investigate  these  effects  on  DM  liquid  water 
accumulation  and  flooding.  The  following  conclusions  have  been 
reached. 

(1)  One-dimensional  model  simulations  demonstrate  that  the 
effect  of  immobile  liquid  saturation  is  critical.  It  significantly 
alters  the  amount  of  liquid  water  accumulation  and  flooding 
level  in  the  DM.  The  effect  is  diminished,  however  when  the 
DM  faces  a  high  liquid  saturation  level,  particularly  with  inter¬ 
facial  liquid  coverage  on  the  DM  surface.  Nevertheless,  the  effect 
of  immobile  liquid  saturation  should  be  considered  in  macro¬ 
scopic  two-phase  models  to  understand  better  the  complex 
water  condensation/evaporation  processes  in  the  DM  and  to 
predict  more  accurately  single-  and  two-phase  regimes  of  the 
DM,  which  are  key  aspects  for  the  water  management  of  PEFCs. 

(2)  When  linear  variation  of  the  cosine  of  the  DM  contact  angle 
is  assumed  along  the  DM  thickness,  the  model  simulations 
showed  that  the  liquid  saturation  profiles  in  the  DM  are  strongly 
influenced  by  the  direction  of  the  DM  contact  angle  variation. 
Liquid  water  transport  is  enhanced  when  liquid  water  is  trans¬ 
ported  from  relatively  hydrophobic  pores  towards  relatively 
hydrophilic  pores.  The  numerical  observations  can  impact  DM 
design  and  optimization  strategies  for  improved  water  man¬ 
agement  of  PEFCs. 

Since  these  two  features  are  of  great  importance  in  DM  liquid 
saturation  profiles  and  flooding  phenomena,  efforts  are  currently 
underway  to  integrate  the  present  one-dimensional  model  into  a 
three-dimensional,  full-cell  PEFC  model  in  order  to  examine  their 
multi-dimensional  effects. 
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